# https://www.jianshu.com/p/9acaa566f04b
#一张图展示多组富集分析结果

library(Seurat)
library(patchwork)
library(clusterProfiler)
# library(org.Mm.eg.db) ##加载小鼠
library(org.Hs.eg.db) ##加载人类
library(tidyverse)
# 导入注释好的pbmc数据集
pbmc <- readRDS(sub_cluster_id, file = "./data/temp/Fibroblast_3亚群_i-m-myCAF.rds")
pbmc$sub_type <- Idents(pbmc)

# 寻找每个群的marker基因
table(Idents(pbmc))
# Idents(pbmc) <- 'sub_type'#因为这个是单样品数据集，在做多样品的时候，把Idents设置为不同样品即可计算不同样品的差异基因
# 
# iCAF  mCAF myCAF 
# 4580   966  3815 

markers <- FindAllMarkers(pbmc)
sig_dge.all <- subset(markers, p_val_adj<0.05&abs(avg_log2FC)>0.15) #所有差异基因
#saveRDS(sig_dge.all, file = "deg.rds")
View(sig_dge.all)

##1. iCAF的上调基因####
sig_dge.iCAF_up <- subset(sig_dge.all, avg_log2FC>0.5&cluster=='iCAF') #iCAF细胞的上调差异基因
#富集分析
ego_iCAF <- enrichGO(gene          = row.names(sig_dge.iCAF_up),
                    #universe     = row.names(dge.celltype),
                    OrgDb         = 'org.Hs.eg.db',
                    keyType       = 'SYMBOL',
                    ont           = "ALL",  #设置为ALL时BP, CC, MF都计算
                    pAdjustMethod = "BH",
                    pvalueCutoff  = 0.05,
                    qvalueCutoff  = 0.05)
ego_iCAF <- data.frame(ego_iCAF)
ego_iCAF$Group="iCAF"

##2. mCAF的上调基因####
sig_dge.mCAF_up <- subset(sig_dge.all, avg_log2FC>0.5&cluster=='mCAF') #mCAF细胞的上调差异基因
#富集分析
ego_mCAF <- enrichGO(gene          = row.names(sig_dge.mCAF_up),
                     #universe     = row.names(dge.celltype),
                     OrgDb         = 'org.Hs.eg.db',
                     keyType       = 'SYMBOL',
                     ont           = "ALL",  #设置为ALL时BP, CC, MF都计算
                     pAdjustMethod = "BH",
                     pvalueCutoff  = 0.05,
                     qvalueCutoff  = 0.05)
ego_mCAF <- data.frame(ego_mCAF)
ego_mCAF$Group="mCAF"

##3. myCAF的上调基因####
sig_dge.myCAF_up <- subset(sig_dge.all, avg_log2FC>0.5&cluster=='myCAF') #myCAF细胞的上调差异基因
#富集分析
ego_myCAF <- enrichGO(gene         = row.names(sig_dge.myCAF_up),
                      #universe     = row.names(dge.celltype),
                      OrgDb         = 'org.Hs.eg.db',
                      keyType       = 'SYMBOL',
                      ont           = "ALL",  #设置为ALL时BP, CC, MF都计算
                      pAdjustMethod = "BH",
                      pvalueCutoff  = 0.05,
                      qvalueCutoff  = 0.05)
ego_myCAF <- data.frame(ego_myCAF)
ego_myCAF$Group="myCAF"

# 分数转化小数
fenshu2xiaoshu<-function(ratio){
  sapply(ratio,function(x) as.numeric(strsplit(x,"/")[[1]][1])/as.numeric(strsplit(x,"/")[[1]][2]))
}

# ego_iCAF
ego_iCAF$GeneRatio1=fenshu2xiaoshu(ego_iCAF$GeneRatio)
ego_iCAF <- ego_iCAF[order(-ego_iCAF$GeneRatio1), ]
write.csv(ego_iCAF,'data/temp/enrichGO_ego_iCAF.csv')
# ego_mCAF
ego_mCAF$GeneRatio1=fenshu2xiaoshu(ego_mCAF$GeneRatio)
ego_mCAF <- ego_mCAF[order(-ego_mCAF$GeneRatio1), ]
write.csv(ego_mCAF,'data/temp/enrichGO_ego_mCAF.csv')
# ego_myCAF
ego_myCAF$GeneRatio1=fenshu2xiaoshu(ego_myCAF$GeneRatio)
ego_myCAF <- ego_myCAF[order(-ego_myCAF$GeneRatio1), ]
write.csv(ego_myCAF,'data/temp/enrichGO_ego_myCAF.csv')

#all
all <- rbind(ego_iCAF[1:15,],
             ego_mCAF[1:15,],
             ego_myCAF[1:15,])

a=all$Description
a_ego_iCAF=ego_iCAF[ego_iCAF$Description%in%a,]
a_ego_mCAF=ego_mCAF[ego_mCAF$Description%in%a,]
a_ego_myCAF=ego_myCAF[ego_myCAF$Description%in%a,]
all <- all <- rbind(a_ego_iCAF,a_ego_mCAF,a_ego_myCAF)
write.csv(all,'data/temp/enrichGO_ego_all_CAF.csv')

# all$GeneRatio <- c(0.3798883,0.3798883,0.3854749,0.2105263,0.1184211,0.09210526,0.1502591,0.1450777,0.0880829)
# 转化以下，不然显示的时候是分数形式
# as.numeric(strsplit("1/8","/")[[1]][1])/as.numeric(strsplit("1/8","/")[[1]][2])
# fenshu2xiaoshu<-function(ratio){
#   sapply(ratio,function(x) as.numeric(strsplit(x,"/")[[1]][1])/as.numeric(strsplit(x,"/")[[1]][2]))
# }
# as.numeric(strsplit(all$GeneRatio,"/")[[1]][1])/as.numeric(strsplit(all$GeneRatio,"/")[[1]][2])
# all$GeneRatio1=fenshu2xiaoshu(all$GeneRatio)


library(forcats)
all$Description <- as.factor(all$Description)
all$Description <- fct_inorder(all$Description)

p1 <- ggplot(all, aes(Group, Description)) +
  geom_point(aes(color=p.adjust, size=GeneRatio1))+theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.x=element_text(angle=90,hjust = 1,vjust=0.5))+
  scale_color_gradient(low='#6699CC',high='#CC3333')+
  labs(x=NULL,y=NULL)+guides(size=guide_legend(order=1))
pdf("./data/output/CAF_GO亚群合并气泡图.pdf",width = 6,height = 10)
p1
dev.off()
